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We use the overlapping grids method to construct a fourth order accurate discretization of a 
first order reduction of the Klein-Gordon scalar field equation on a boosted spinning black hole 
blackground in axisymmetry. This method allows us to use a spherical outer boundary and excise 
the singularity from the domain with a spheroidal inner boundary which is moving with respect to 
the main grid. We discuss the use of higher order accurate energy conserving schemes to handle the 
axis of symmetry and compare it with a simpler technique based on regularity conditions. We also 
compare the single grid long term stability property of this formulation of the wave equation with 
that of a different first order reduction. 



I. INTRODUCTION 

Black hole excision has become an important tech- 
nique in numerical relativity. First proposed by Unruh 
P|, excision consists in placing an outflow inner bound- 
ary which eliminates the black hole singularity from the 
domain. Sometimes combined with singularity-avoiding 
slicings, it is used in all of the current long-term black 
hole evolutions. (For review articles see Refs. 
some more recent work includes 0,0,I3-) Excision has 
proved most successful with black holes at a fixed co- 
ordinate location. For example, in studies of orbiting 
compact objects (black holes and/or neutron stars), long 
runs have been achieved in co-rotating frames, where dy- 
namic gauge conditions attempt to keep the black holes 
at fixed locations. For more general orbits, and for sim- 
ulations over many orbital time scales, it is anticipated 
that the ability to move the black holes on the compu- 
tational grid through a stable excision algorithm will be 
crucial. 

Excision methods for moving black holes typically re- 
quire the extrapolation of data onto the trailing edge of 
the black hole 

B 113, in El El- We recently reported 
on a new excision method using simultaneous, multiple 
coordinate patches, and implemented these numerically 
using overlapping grids [T3 | . In this method each bound- 
ary, whether an outer boundary or an excision bound- 
ary, moving or static with respect to a main coordinate 
system, is fixed to at least one coordinate system. The 
different coordinate patches overlap just like the charts 
of an atlas. Information is communicated from one grid 
to the other using only interpolation, without any de- 
composition into ingoing and outgoing variables. This 
excision method has a number of advantages, including: 
(1) the possibility of choosing coordinate patches which 
conform to each individual boundary, giving smooth nu- 
merical boundaries and thus simplifying the boundary 
treatment; (2) the ability to adapt the coordinates near 
the black hole to the horizon geometry, allowing the ex- 
cision surface to be placed relatively far from the sin- 
gularity, and thereby "excising the excisable"; (3) the 
use of simple data structures and the implementation of 



standard methods for parallelization and adaptive mesh 
refinement due to the fact that individual grids are logi- 
cally Cartesian. 

In our previous work we demonstrated this excision 
method by evolving a Klein-Gordon massless scalar 
field on a boosted Schwarzschild background. Numer- 
ical tests showed that the scheme was stable and sec- 
ond order convergent for very large boost parameters 
(v = 0.98c). Long-term convergence test of spherical 
waves in Minkowski space also showed that the interfaces 
between the overlapping grids did not introduce growing 
numerical errors. Moreover, our tests indicated that the 
overlapping grids technique is robust, in the sense that 
stability was not dependent on some of the fine details of 
the implementation. For example, stability was indepen- 
dent of the interpolation method, the number of grids 
used, the physical overlap size of the grids (which was 
kept fixed while testing convergence), and their relative 
resolutions. 

Overlapping grids have also been used recently in evo- 
lutions of the full Einstein equations. Thornburg for 
example, used six spherical coordinate patches to excise a 
Kerr black hole (a = 0.6). The patches arc designed such 
that coordinate lines in one direction overlay each other, 
so that the interpolation required is effectively one di- 
mensional. Using the BSSN formulation, Thornburg was 
able to evolve the Einstein equations for 1500 Af , and it 
appeared that the instabilities were related to the outer 
boundary treatment. He also used overlapping patches 
in an apparent horizon solver, which is now available in 
Cactus [il. 

More recently, Anderson and Matzner ^3 have per- 
formed simulations of black holes with overlapping grids. 
They used two stereographic patches to cover the spher- 
ical excision region, and solved the standard g-K ADM 
equations in a constrained evolution. They were able to 
evolve boosted black holes [v ~ 0.5c) across the com- 
putational domain, and runs of single, stationary black 
holes ran for more than 700 M . 

Finall y, in a related approach, Reula, Tiglio and 
Lehner [l^ use multiple coordinate patches designed 
such that all boundary points on neighboring patches 
arc aligned. Rather than overlapping, the grids just 
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touch. The communication between the touching grids 
is based on the Carpenter-GottUeb-Abarbanel's method 
which requires the computation of the characteris- 
tic variables. They successfully evolved a Klein-Gordon 
scalar field on a Kerr background spacetime and are cur- 
rently working towards fully relativistic evolutions. 

In this paper we expand upon our previous work, and 
choose to study improvements in the following directions: 
(1) we generalize the background spacetime geometry to 
allow for spinning black holes, and investigate the effect 
of high spin (a = 0.99) on the system's stability; (2) 
we upgrade our derivative and interpolation operators 
to give global fourth order accuracy; (3) we introduce 
a spherical patch for the outer boundary, removing all 
boundary corners. 

While in we put strong emphasis on the use of 
energy conserving schemes based on difference operators 
satisfying the summation by parts rule, here we allow 
ourselves to explore discretizations for which a stability 
proof based on the discrete energy method is not imme- 
diately available. At times, such discretizations can be 
much simpler, particularly in the higher order accurate 
case. 

The paper is organized as follows: In Sec. ^ we out- 
line the main ingredients of the overlappin g gr ids method. 
Scc. lIIIl is a generalization of Sec. III.B of to the spin- 
ning case with the inclusion of a potential in the wave 
equation. The various coordinate systems used and the 
regularization of the equations on the axis are discussed 
in Sec. IIVI Before describing in Sec. IVII the details of 
the discretization, in Sec. we produce group veloc- 
ity diagrams which illustrate how information propagates 
through the domain. A high resolution convergence test 
using the forcing solution method is described in Sec. lVIII 
This paper contains three appendices. 



II. OVERLAPPING GRIDS 

When solving hyperbolic initial-boundary value prob- 
lems numerically, it is often difficult, if not impossible, 
to accurately represent the entire domain, boundary in- 
cluded, with a single grid. For example, at least two 
coordinate patches are required to cover the entire sur- 
face of a 3-sphere without coordinate singularities. When 
solving for the spacetime representing a binary black hole 
collision, one may want to use a spherical domain with 
the outer boundary placed sufhciently far away, in which 
the black hole singularities are excised. To represent the 
interior of a sphere with two (or one after the merger) 
smaller spheres removed, more than one coordinate sys- 
tem is needed, specially if one wants the coordinates to 
be adapted to the boundaries of the problem. 

The overlapping grids methods provides a simple and 
flexible solution to such problems. Our scheme is based 
on that described in Starius [23| and Ref . . To demon- 
strate how the algorithm works we consider the first order 




FIG. 1: This illustrative example consists of two overlapping 
grids. The one on the left is the main grid. The other one 
is adapted to the physical boundary of the problem, which 
may be moving with respect to the main grid. The gridpoints 
marked with a solid dot are updated via interpolation. The 
overlapping grids algorithm is described in the body of the 
paper. 

linear hyperbolic system 

dfU ^ P{t,x,dj:)u , (I) 

where m is a vector valued function representing tensor 
field components, in a domain Q{t) with smooth bound- 
ary 9f2(t). The problem includes appropriate initial and 
boundary data. In general, system JQ) is specified in 
a coordinate system {t, x} which is not and cannot be 
adapted to the boundary of the problem. We will re- 
fer to this boundary as the physical boundary, and we 
allow for it to be moving with respect to the main coor- 
dinate system. We introduce a coordinate system {t', x'} 
adapted to d^l{t), in the sense that the boundary sur- 
face can be represented by x'* = const., for some i. The 
equations in this coordinate system will take a different 
form 

dt'u' = P'{t',x',d,,)u' , (2) 
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where the components of u' are related to the compo- 
nents of u via tensor transformation laws. The initial 
and boundary data are also transformed. To maintain 
simultaneity of patches we restrict the coordinate trans- 
formation by demanding that t' = t. 

The problem is discretized by introducing a grid for 
each coordinate system. Any grid point of the main grid 
which lies beyond a x'* — const, line (or surface) are 
dropped, as illustrated in Fig. ^ The two grids over- 
lap and the physical overlap size is kept constant when 
performing grid refinement tests. 

The numerical computation of the right hand side of 
at a particular gridpoint Xij requires information 
from a number of gridpoints in each coordinate direc- 
tion. With a fourth order accurate centered difference 
operator, for example, one needs two gridpoints in each 
direction. Gridpoints of the main grid, at which the right 
hand side cannot be evaluated, are updated via interpo- 
lation from the other grid. All components of u' are in- 
terpolated onto the main grid and the tensor law trans- 
formation is used to evaluate the components of u. A 
similar procedure is done at the non-physical boundary 
of the second grid. Points that require interpolation are 
marked with a solid dot in the figure. The numerical 
treatment of the physical boundary is as that of a single 
grid boundary. 

In our code we use n-th order Lagrange interpolation, 
which in two dimensions is given by 

n n n n 

i=l 7 = 1 1 = 1 ' ' k=l "J "'^ 

If / is sufficiently smooth, the interpolating function is a 
n-th order approximation of f{x,y). 

Note that if the physical boundary moves with respect 
to the main grid, gridpoints may have to be dropped from 
or added to the main grid and the set of points which re- 
quire interpolation needs continuous updating. This is 
done at the end of each full time step of the time integra- 
tor of choice (e.g. fourth order Runge-Kutta). When grid 
points are added to the main grid, the grid adapted to 
the boundary is able to provide accurate data for these 
points. 

The overlapping grid s method requires artificial dissi- 
pation for stability [221 . We use sixth order dissipation, 
see Eq. (|52f) . which has a seven point stencil in each direc- 
tion. This means that in our code we actually interpolate 
three, rather than two gridpoints. 

Here we consider the case in which a scalar field propa- 
gates on a boosted spinning black hole background. Two 
boundaries are introduced: an inner and an outer bound- 
ary. The first one represents the excision surface and 
is purely outflow. It requires no boundary data. The 
second one is introduced for computational reasons. We 
need to have a bounded spatial domain because of limited 
computational resources. To handle the two boundaries 
we introduce two additional coordinate patches. One 
patch is adapted to the outer boundary and one patch is 




FIG. 2: We use a main cylindrical grid and two spherical 
grids adapted to the inner and outer boundaries. The irreg- 
ular shape of the main grid is the result of having dropped 
gridpoints which lie beyond imaginary lines on the spherical 
grids. The dot represents the ring singularity of the black 
hole. 



co-moving with the black hole, and fixed to the inner exci- 
sion boundary. We choose cylindrical coordinates for the 
main coordinate patch. The black hole is boosted with 
velocity /3 along the axis of symmetry with respect to 
this coordinate system. Spheroidal coordinates are used 
on the second patch, such that the location of the event 
horizon is at a constant coordinate value, and spheri- 
cal coordinates on the outer patch. We require that all 
data in all three coordinate systems be simultaneous. By 
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adapting these coordinates to the black hole horizon, we 
may excise the spherical grid at the event horizon for all 
values of the boost parameter. This allows for an effi- 
cient use of the excision technique, as we can excise the 
excisable. The layout of the grids is illustrated in Fig. |5] 
For the fourth order accurate case, which is the minimum 
order that we require, Eq. (O on an uniform grid takes 
the form 



/int(a;i + aAx, yj + bAy) = 

1+2 J+2 i+2 j+2 

E-^ -pr a + i- l -pr 6_K7__fc 
2^ 11 n-7 11 - ■'P' 



(4) 



q=j 



-1 l = i~l ^ k = j-l ^ 



where < a < 1 and < < 1 . 

The boundary treatment and the discretization near 
the axis of symmetry are described in Sec. IVII 



The characteristic speeds in an arbitrary direction n, 
with |n| = 1, are given by the eigenvalues of the matrix 



^" = A'n, 












(10) 



-7") 



These are s± ~ (7*" 

/3" ± aVh/^ and sq — with multiplicity equal to the 
spatial dimension of the problem, where a is the lapse 
function, /?' the shift vector, and hij is the induced 3- 
metric on the t = const, slices in the Arnowitt-Deser- 
Misner (ADM) decomposition ljCl|l . Hyperbolicity re- 
quires that the characteristic speeds be real, namely, 
(^t")2 _ ^tt^nn ^ det(/i.y)ft."" > for any n, which 
will be true as long as the t ~ const, hypersurfaces are 
spacelike. 

One can verify that 



III. A FIRST ORDER REDUCTION OF THE 
WAVE EQUATION 

In this section we write down the wave equation on a 
general curved background and discuss a particular first 
order reduction. In Appendix lUl we consider a different 
first order reduction, which has different features, and 
compare the numerical stability of the two formulations. 
For a definition of strong and symmetric hyperbolicity 
we refer the reader to |2ll [23 . The rest of this section is 
a generalization of Sec. III.B of [T^ . 

The equation of motion for a scalar field propagating 
on a curved background (M, g) is given by the second 
order wave equation 



(5) 



where V denotes the covariant derivative associated with 
the Lorentz metric g and V{^) is a potential. 

In terms of the tensor density 7^'' — where 
g = det{gfj,^), the wave equation can be written as 

If wc introduce the auxiliary variables T = dt^ and di 
di^, we can rewrite Eq. © as a first order system, 

dt^ = T, 

dtT = -(7"a,r + 9,(7''*T) + c),(7*^rf,) 

, dV 



did. 



(6) 

k = 

(7) 
(8) 

(9) 



An attractive feature of this particular first order reduc- 
tion is that the constraint variables propagate trivially, 
namely df Ci — 0. This ensures that any solution of 
which satisfies the constraints initially, will satisfy them 
at later times, even in the presence of (static) boundaries. 



Hit,x)^ I -777 




^00 

tt 

?77 



(11) 



where ^ and rj are functions of t and x, is the most 
general symmetric matrix that satisfies HA^ = {HA^)'^ , 
i = 1, 2, 3. When positive definite, which will be the case 
if and only if dt is timelike and ^ > 0, 77 > 0, it represents 
the most general symmetrizer of system (|7| IHl . Using the 
fact that 7** < (because the t = const, slices are space- 
like) and 7*^ ^ V^9'^ = ay^M(h~){h'^ - 
one can show that the symmetrizer is positive definite if 
and only if the vector field dt is timelike. 

The symmetrizer can be used to construct an energy 
and obtain energy estimates. The time derivative of 



E 



{-j''T^+f^d,dj) + V^V{<^) 



is given by 
dt 



(Tj^'T + Tf^dj) m d^a 



d^x (12) 



(13) 



{Tda^'T + 2Tda'^dj + dAY^dj) d^x 



I dtV^V{'P)d^x, 



where rii is the outward pointing unit normal to dfl. If 
y($) is quadratic in <!>, e.g. V — im^$^, dt is time-like, 
and maximally dissipative boundary conditions are used, 
we have an energy estimate. Note that in this case (|12|) 
corresponds to the choice ^ — m?/2 and = 1/2 in Hll|l . 
If, furthermore, the background admits a time-like vector 
field k and wc use a coordinate system adapted to it, 9t = 
k, the components of 7'''' will be time independent and 
we have, ignoring boundary terms, energy conservation. 
The integrand of the surface term can be written as 



2iTj*'T+Tj'^dj)ni 



(14) 
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where A± = 7" ± 7*" and 



and r is determined implicitly by 



,(±A±;«) 



7™d, 



^/2 



\/2 v/T^ 



(0;n) 



(15) 
(16) 



2 2 



or explicitly by 



7 



are the orthonormal characteristic variables of HA^ . To 
simplify the notation we have introduced the quantities 
yS~7^™7^, 7''" = 7'^"/7" and 7^. The latter 

«i7l7i ~ 1 ^'^'^ = 0. To express the 

primitive variables in terms of the characteristic variables 
we invert Eqs. ((T3|l - lfTC|) . 



' BL 



(23) 

where = + + z^. The inverse metric can be 
written as 



9'^ 



(24) 



T 



Vl+7*" ...(+A.;„) _ Vl-7*" ...(-A_;„) 



^/2 



^/2 

y^(-A_;n) 



(17) 



\/2 wi + 7*" 



+ 7lti;(°'")(18) 



We use Eqs. (|15|) -(I18 |I in the boundary conditions to pre- 
scribe exact data to the incoming characteristic variable 
The variable <I> is a zero speed characteristic 



w 



(+A+;n) 



variable for any direction n and requires no boundary 
data. 

We also assume axisymmetry, which implies that there 



where = rj'^'^ip — g^'^lv is a null vector. 

The quantities M and a are constants, M representing 
the mass and Ma the angular momentum of the black 
hole as measured from infinity. We restrict ourselves to 
the case a? < AP. We recall that the event horizon 
is located at r = r+ = M + {AP - a^)^/^, the Cauchy 
horizon at r = r_ = AI —{M^ — a^y^^ , and the stationary 
limit surface, the set of points in which the Killing field 
K = d/dt becomes null, is given by r = M + (A/^ — 

cos^ 6*)^/^. Another set of points in which k is null is 
given by r = M- (M^ - cos^ 61) ^Z^. The disc a:^ + 



< 



corresponds to r = 0. The ring x 



y 



exists a spacelike Killing field rp = ip'^df, = d^. By adopt- introduce the quantity p|l(7-. 



is a curvature singularity. For later convenience we 



ing coordinate systems adapted to the Killing field, we 
have that the metric components are independent of the 
(j) coordinate. Since we arc only interested in axisymmet- 
ric solutions of the wave equation, i.e., solutions which 
do not depend on the variable d^, which represents 
9^$, can be eliminated from the system. 

In the next section we define the various coordinate 
systems used in our overlapping grids scheme and spe- 
cialize Eqs. (0-10 to these coordinates. 



IV. THE KERR METRIC IN KERR-SCHILD 
COORDINATES 

We are interested in the case in which the background 
is given by a rotating Kerr black hole. We will use Kerr- 
Schild coordinates [24 . |25| . After recalling basic proper- 
ties of these coordinates we explicitly compute the coef- 
ficients 7'^'^ needed in the evolution system and 
determine the regularized equations on the axis. 

The Kerr-Schild metric components are given by 



A. Boosted cylindrical coordinates 

The main coordinate system is obtained by performing 
a Lorentz boost, followed by a transformation to cylin- 
drical coordinates. Under a Lorentz boost, i.e., in the 
new coordinates 



t 

X 



lit - dz) 

X 



(25) 



y ^ y 

z = j{z-pt), 



where 7= (1-/3^) the components of the Kerr- 
Schild metric become 



77^p + 2ff^p4 , 



jypp = diag{-l,+l,+l,+l}, 
' f rx + ay ry ^ ax z 



where 



77^^ = diag{-l,+l,+l,+l} , 
H 



2^2 ' 



1, 



rx + ay ry 



ax z 
r 



a 



(19) 

(20) 
(21) 

(22) 



where f — 7(r 4- /3z), z ~ 7(2 + /3r), z ~ 7(2 -|~ (3t) and 
r = rBhix,y,z) = rBLix,y,iiz + pt)). At time i the 
singularity is located at x^ + y^ = a^, z = —j3t. 

Wc_now go to cyhridrical coordinates {t,p,z,<f)}, with 
pcos0 = X and psmcj) = y, obtaining 



rp 



ap 



r r'^ + a^ 



+ a^ 



diag{-l,+l,+l,+p2} 
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and 



diag{-l,+l,+l,+— }^ 
P 



rp 



+ a^^ r' 



(26) 



where z — ^{z + pi) and r ~ tbl (p cos 0,p sine/), 7(2 + 
Pt)). Now the singularity is located at p = |o|, z = — /3f. 
Note that in these coordinates the time derivative of ^/— g 
vanishes. 

In this coordinate system the symmetrizer associated 
with this first order reduction does not lead to a con- 
served energy (dt^^'^ 7^ 0), except for P ~ 0. The region 
in which the system is symmetric hyperbolic is deter- 
mined by the set of points in which dt is timelike, namely 



9tt 



1 ^ > 0. 



(27) 



In Sec. El we analyze the group velocity of the system 
and produce plots which give a graphical representation 
of how information propagates throughout the domain. 



B. Co- moving spherical coordinate system 

To excise the black hole we introduce a spherical co- 
ordinate system {t' ,r' ,9' , (f>'} which is related to the 
unboosted Cartesian coordinates {t,x,y,z}, the coordi- 
nates in which the black hole is at rest, by 



t' = 



i^j{t-pz) 
rB-L{x,y,z 



(28) 



~ tan 



tan 




cos 



\rBi,x + ay 



In the unboosted case the system is symmetric hyperbolic 
outside the stationary limit surface, and only strongly 
hyperbolic inside (see Fig. 0J. For P ^ the region of 
lack of symmetric hyperbolicity increases. This issue was 
explored in greater detail in [Tj| for the a = case. 

On the axis of symmetry (p = 0) equation (jS)) needs 
to be regularized. We want to express it in a form which 
avoids "0/0". This can be done by taking the limit p — > 
in the equations. For this purpose it is convenient to 
introduce the rescaled quantities 







-9 ' 


P 




P^ 


ryPP 




ryP^ 


P 







which have a finite limit for p provided that r > 
To > 0. The right hand side of (|SJ) at p = becomes 



Note that for /? = the time coordinate coincides with 
the Kerr-Schild time and not that of the Boyer-Lindquist 
coordinates [2^. In addition, the azimuthal angle 0' does 
not coincide with the Boyer-Lindquist (jSibl, but it is re- 
lated to it via 



+ a / (r^ - 2Mr + a^y^dr 



by 



The inverse coordinate transformation of H28() is given 



t = 

X = 

y = 

z = 



Pr' cos e' 



7 

sin 0' (r' cos </>' — a sin 0') 
sin 6' (r' sin (p' + a cos 0') 
r' cos 9' 



(29) 



The coordinates are adapted to the event horizon in the 
sense that its location (r' = r+) is time independent and 
setting t' =t ensures simultaneity with the main coordi- 
nate system. The components of the inverse metric can 
be constructed from 



= (-7(1 + /? cos 6*'), 1,0,0) 

/ -Pbl -Pl{r'^ + 0."^) cos 0' p-fr'ainB' -a/37 cos 6*' \ 



nP " 



v ■ 



-2 a' 



(30) 
(31) 



where p^^ 



r'^ + a? cos^ 9' . By definition, 7'^ " 



, where 5^ = if " - 2H^^' t and 
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Note that for a 7^ there is a region outside the event 
horizon in which the system is not symmetric hyperbohc, 
even when f3 ~ 0. 

As we did in the cyhndrical case on the axis of symme- 
try {9' = or 9' = n) we need to regularize the equations. 
Introducing the quantities 



-ft 

7 



7 



t't' 



sm 
sin 6 



7, f'^' 



7 



t'r' 



sm 
e'9' 



7, f'^' 



7 



t'e' 



sin^ 9' 



and taking the hmit 9' 6^, where 9q — Q,tt, the right 
hand side of ^ on the axis becomes 

dt,T' = (f '-'a^-T' + a^'Cf '''T')±2f' ^'t' (32) 

+dr'{r''''dr') + 2j'>'''dg,dg, (33) 

where the components of 7^* " are understood to be eval- 
uated at 9' = 0,TT. 

Note that in general the characteristic speeds in the 
rotating case are higher than those of the non-rotating 
case. For example, in the radial direction at (r, 9) = 
(M, 7r/2) they can be as large as as 4/3 in the extremal 
unboosted case a = M, /3 = 0. In the nonrotating case, 
instead, the same characteristic speeds are bounded by 1 
outside the event horizon. 



C. Outer spherical coordinate system 

In order to have a smooth outer boundary we introduce 
a spherical coordinate system {t, f , 9, </>}, which is related 
to the boosted cylindrical coordinate system through the 
time independent transformation 



P 

z 



r smt 
f cos ( 



and its inverse 



= tan"^ - 



which ensure simultaneity. 

In this coordinate system we have 



(34) 



(35) 



As for the co-moving spherical coordinate system, to 
have the evolution equation in a form that avoids "0/0" 
it is convenient to introduce the rescaled quantities 







rytr 






sin^' 


sin 9 ' 


sin^e' 




f 








sin ^' 


sin^ 9 


sin 9 



We get 



dif = (^f'^dff + df{Y'T)±2f^f 

+a,-7"T + d-^rd, - /(-7") • 



D. Transformation of variables 

The interpolation procedure which is used to commu- 
nicate information between the grids requires the co- 
ordinate transformations and the tensor transformation 
laws to relate the components of the evolved fields. The 
boosted cylindrical and co-moving spherical coordinate 
systems are related by 



t = t' 

p = s\n9' \/ 7''^ + a? 
z = -f-^r' cos 9' - (3t' 



(36) 



tan 



and the inverse transformation 



t' = t (37) 
r' = rBLipcoscj), psin(j),j{z + f3t)) 

9' = tan-i ^ ^"^^ ^ 



z 



= COS ^ y 

h' = 0-tan"^ 



z 

V^BL 



a 



where z = 7(2 -I- (it). 

The field $ transforms like a scalar and the fields T and 
di transform like components of 1-forms. In this case we 
have 



-5s 



sin ( 



diag<^ 



rr 



r r'' + a'' 



■ sm 



sin^ I 



z 

— cost 
r 



— - sm cos U sm tl, — 

rr + 



dp 



T' + -f/3- 



■ cos 9'dr' -jP—sin9'dg{38) 



r'\/r'^ + a 



r'2 -t- 



Pbl 

sin 9'dri H ^ cos 9' do' 



pIl 



7 5 cos 9'dr' — 7^5— sin 9'de' , 

Pbl Pbl 
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and 



dr' 



(39) 



sin 6' dp + J ^ cos 9' dz 



dg, = \Jr''^ + a2 cos 9' dp - 7 Vsiiifl'd^, 

where {T,dp,dz) and {T',dr',dei) are the fields on the 
cyhndrical and spherical grids, respectively. 

We also need the transformations between the boosted 
cylindrical and the outer spherical coordinate systems, 
Eqs. (El 123), and 

— COS 

dp = sin —dg (40) 

r 

sin 9 

dz — cos9df —dg (41) 



and 



df = sin9dp + cos 9dz (42) 
dg = f (cos^dp — sin^dj) (43) 



V. GROUP VELOCITY DIAGRAMS 

To have a graphical representation of how information 
propagates through the domain it is useful to construct 
group velocity or light cones diagrams. These diagrams 
are obtained by looking at the group velocity of the sys- 
tem at selected points of the domain. We briefly describe 
how the procedure works. More details can be found in 

Consider a system of quasilinear partial differential 
equations in two spatial dimensions 



dfU = A*(i, X, u)diU + B{t, X, u) . 



(44) 



The characteristic speeds, or phase velocities, in the di- 
rection n = (ni, 712) are given by the eigenvalues of A^rii. 
We focus on one of them, which we denote hy X{t, x,u,n). 
The group velocity is given by its gradient with respect 
to n, i.e., 



/ ^ 1 dX dX 

Vg(t,x,u,n) = - — , -— 
' oni on2 



(45) 



In order to produce light cones plots at a certain time 
^0, about a solution uq, we select a number of uniformly 
spaced points of the domain, Xi, and plot the parametric 
function 



x{t) = Xi + aVg{to,x^,uo,n{T)) 



(46) 



where nlr) = (cost, sin t) with r G [0, 27r) and a is a 
constant introduced for aesthetic reasons. It avoids that 
the cones either overlap or are too small. 

We produced plots for the linear system {TJ)-© in the 
cylindrical coordinate system. Fig. O illustrates the light 



cone structure in the non-rotating case, whereas Fig. 0] 
represents the light cones in the rapidly spinning case, 
and Fig. |31 shows the rapidly spinning, highly boosted 
case. It is interesting to see how in Fig.|Slthe light cones 
tilt. While the light cones that are being approached 
by the black hole are almost unperturbed, those behind 
are tilted in such a way that the information is forced to 
follow the black hole. This "dragging" effect is clearly 
noticeable in our numerical simulations. 
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FIG. 3: The light cone structure in Kerr-Schild cylindrical 
coordinates for P = 0, a = 0, and M = 1. The event horizon 
i2 = 2M. 



is located at r = yp^+T 



VI. DISCRETIZATION 

To discretize system ((TJ-© we replace the partial 
derivative di with the finite difference operator Di , with- 
out expanding derivatives of products. This leads to the 
semi-discrete system 

5t$ = T (47) 
dtT = -{j''D,T + D,{Y'T) + D,{Y^d,)+ (48) 
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FIG. 4: The light cone structure in Kerr-Schild cylindrical 
coordinates for /3 = 0, a = 0.99M, and M = 1. The red line 
represents the event horizon and the blue line the Cauchy 
horizon. The brown dot is the ring singularity. The outer 
region between the event horizon and the green line, the sta- 
tionary limit surface, represents the ergoregion. Between the 
two green lines the vector field dt is spacelike and, therefore, 
the system is only strongly hyperbolic. 



FIG. 5: The light cone structure in Kerr-Schild cylindrical 
coordinates for /3 = -0.85, a = 0.99M, and M = 1. The red 
line represents the event horizon and the blue line the Cauchy 
horizon. The brown dot is the ring singularity, which is mov- 
ing towards positive values of z. The outer region between the 
event horizon and the green line represents the ergoregion. 



dtd, = D:T. 



(49) 



We now discuss the discretization near the axis of sym- 
metry and near the physical boundaries. 



where we have not explicitly written out the gridfunction 
indeces. As was shown in [ij], this type of discretization 
leads to discrete energy conservation provided that dt is 
a time- like Killing field • 

In the interior we use the centered fourth order accu- 
rate finite difference operators 



(1) 



- n(2) 



(1)^(1) 



8u, 



(50) 

-I- Ui_2j)/(12/ll) 

D'^' I u,j (51) 

(-■"i.j+2 + 8Uij + i - 8Uij_i + liij_2)/(12/l2) 



h- „(2) „(2) 



A. Axis of symmetry 

The difference operators (|50f) and (|5I|I have a 5 point 
stencil. Gridpoints that arc too close to the axis require 
special treatment. If a gridpoint is close to but not on 
the axis of symmetry we use Eqs. H47|) - (|49|l combined 
with the regularity conditions to evaluate difference op- 
erators in the direction normal to the axis. With respect 
to the coordinate normal to the axis <&, T, and cLa are 
even and dn is odd, where n is normal to the axis and A is 
tangent. If a gridpoint lies on the axis we use the regular- 
ized equations combined with the regularity conditions. 
Specifically, in the boosted cylindrical coordinate case we 
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use 

In the co-moving spherical grid on the axis of symme- 
try (9 = and 9 = n) wc use the approximation 

7 d^J'^ ^ 
whereas in the outer spherical grid we use 

+an"T + ditd, - --'^)/(-7") 

in) 

The operator Dreg represents the centered fourth order 
accurate operator computed using the regularity condi- 
tions. Note that this discretization is not energy conserv- 
ing. This is discussed further in Appendix ^ Clearly, 
the derivative operator in the direction parallel to the 
axis needs to be modified near the physical boundary 
and boundary data needs to be provided at the outer 
boundary. This is discussed in the next subsection. 

B. Boundary conditions 

We always place the excision surface r' = r'^-^^ between 
between the Cauchy and the event horizon, in which case 
no boundary conditions are required. The fourth order 
accurate difference operator near the excision surface in 
the normal direction to it, however, needs to be modi- 
fied. The modification is given in Appendix IXI A similar 
modification of the difference operator is required near 
the outer boundary of the outer spherical grid, f = fmax- 
In addition, here we give data to the incoming character- 
istic variables in the normal direction, including at those 
points which lie on the axis, where the normal is cho- 
sen to be parallel to the axis (see Fig. 4 of dl). We 
do not couple the ingoing to the outgoing characteris- 
tic variables. Our numerical implementation is based on 
Olsson's method |30l |. 

C. Artificial dissipation 

It is known that ove rlap ping grids require artificial dis- 
sipation for stability |23|. Therefore to the right hand 



side of the discretized system we add a term of the form 

Qdu.^ = o (h\{D^l^D^-^f + \A{DfD^^^f) u,, . (52) 

We modify this operator near the physical boundaries 
according to 

QdUQ = ah^D^UQ , 

QdU2 = cr/i^ [Dl - iD^D\ + iD+D'i) U2 , 

where we have only kept the relevant gridpoint index. 
The modification near i = is obtained from the above 
with the replacements i ^ N — i and D± —Dzp. This 
modification was derived by combining the extrapolation 
conditions h^D^u^j = and h^D^UN+j, J = 1, 2, 3 with 
(I52|l whenever gridpoints outside the domain are needed. 
The result is then multiplied by h to ensure that the 
modification is a third order correction, which does not 
lower the global accuracy of the scheme [sj • Near and 
on the axis of symmetry dissipation (in the normal di- 
rection) is computed exploiting the regularity conditions 
of the fields. We were not able to modify the dissipative 
operator such that a discrete energy estimate holds and 
without lowering the global accuracy of the scheme. 

D. The discrete energy method 

Unlike for the second order accurate case considered 
in we are not able to show that our discretization 
satisfies a discrete energy estimate, not even on a single 
grid and with homogeneous boundary data. We should 
stress that the inability to obtain a discrete energy esti- 
mate is by no means a proof of instability. The discrete 
energy method is only a sufficient condition for stability. 

In Sec. IVIII numerical experimentation is used to de- 
termine the stability of our scheme. 

E. Choice of time step and dissipation parameter 

We obtain the fully discrete system by integrating the 
semi-discrete system of ordinary differential equations 
with fourth order Runge-Kutta. Whenever explicit fi- 
nite difference schemes are used to approximate hyper- 
bolic problems, the ratio between the time step size k and 
the mesh size h = imn{hi}, the Courant factor, cannot 
be greater than a certain value j32l |. which is inversely 
proportional to the characteristic speeds of the system. 
In Fig. we numerically estimate allowable values for 
the Courant factor by examining the 2D wave equation 
written in first order form, dtUo = d^Ui, dtUi = diUg. We 
plot the Courant limits for fourth order Runge-Kutta as a 
function of the artificial dissipation parameter, assuming 
fourth order, centered differencing for the spatial deriva- 
tives (|50f) - (l51|l and sixth order dissipation H52|) . A simi- 
lar plot for the second order accurate case can be found 
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in Fig. 7 of [1J|. From the plot we see that the value 
a ~ 0.0075, while damping the high frequency modes, 
docs not force us to reduce the Courant factor. This 
is the value of dissipation parameter that we use in our 
simulations. It is about three times smaller than the 
correspondent value for a second order accurate approx- 
imation (cr w 0.02). 



^ 1.0 




FIG. 6: The maximum value of the Courant factor for the 
stability of the fourth order accurate approximation of the 
2D wave equation, dtuo = diUi, dtUi = diUo, with artificial 
dissipation 15211 . integrated with fourth order Runge-Kutta, 
as a function of the dissipation parameter a. 



VII. NUMERICAL EXPERIMENTS 

In this section we outline the test carried out to verify 
that our overlapping grid code is fourth order convergent. 
We faced some unexpected issues when testing for long 
term stability. These are discussed in Appendix lUl 

Let u be the exact solution of the continuum problem 
at time t and v'^ the solution of the fully discrete approx- 
imation at time step n, such that t ~ kn, obtained with 
a mesh size h. The convergence rate is computed as 



hi = log 



\\u-v''\\h 



(53) 



where || • \\h is a discrete L2-norm. Neglecting round-off 
errors one should expect that lim/i^o Qh = 4 for a fourth 
order accurate scheme. To use this equation wc need an 
exact solution of the continuum problem. We use the 
same forcing solution technique described in [l^ . 

Let us rewrite the partial differential equation as 
L{u) ~ and let w be an arbitrary function. If w is 
inserted into the equation, in general, it will produce a 
non vanishing right hand side. 



Clearly, the modified equation L{u) = L{u) — F = has 
w as an exact solution and the convergence of the code 
can be tested using Eq. (|^ . 

We chose w{t, r, 6) = sin{t + r) cos{n9), where {t, r, 6} 
are the spheroidal coordinates of the rest frame and n is 
an integer. This is an exact solution of 



aw 



(55) 



as long as F is given by 
F = 



sin^ I 



■ sin(t + r) cos{n9) (56) 



2r 

-^5— cos(t + r) cos(n9) 
Pbl 
„ a 

■ sin(t + r) sin(n6 



_dV 
dw ' 

Both w and F are scalar quantities. The evolution equa- 
tion © is modified according to 

dtT ^ -(^^''d,T + d,{f'T)+d,{Y^dj) (57) 

, dV 



L{w) = F. 



(54) 



where g = det((7^j/). On the axis of symmetry we use the 
limits lime^o^i^ = n and lime-,^ = (-l^+'^n. 

The result of our convergence test for M = 1, /? ~ —0.85, 
a ~ 0.99, and n = 2, confirms that wc have fourth order 
convergence on each grid and is shown in Fig. \7\ The 
simulation was stopped when the inner spherical grid was 
about to touch the outer spherical grid. 

Furthermore, Fig. |S1 illustrates that failing to add ar- 
tificial dissipation to the scheme can lead to numerical 
instabilities. 



VIII. CONCLUSION 

The difficulties associated with the construction of a 
stable scheme for the wave equation on a black hole back- 
ground are also present in simulations of fully nonlin- 
ear general relativistic systems. Therefore, the ability to 
handle a scalar field propagating around a boosted spin- 
ning black hole is one of the first necessary steps towards 
the construction of a successful, long term stable binary 
black hole collision numerical code from which gravita- 
tional wave templates can be extracted. In this work 
we showed that the overlapping grids method is com- 
patible with fourth order accuracy, provided that sufla- 
cicntly high order interpolation and artificial dissipation 
are used. We also noticed that simple schemes at times 
work better than more complicated techniques based on 
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conservation of discrete quantities (see, for example, the 
treatment of the axis in Appendix^. 

In Appendix |0 we observed that different hyperbolic 
first order reductions, which have different levels of hy- 
perbolicity (symmetric or only strongly hyperbolic) in- 
side (or near) the excision region, have different stability 
properties. This was unexpected since in the continuum 
limit the treatment of the excision region, if contained in- 
side the black hole, should not affect the solution outside 
and therefore should not introduce any instabilities [48l|. 
We noticed that for the T- formulation the behavior of the 
solution is highly sensitive to the discretization (chang- 
ing the way the dissipative operator is modified at the 
inner boundary, for example, can have noticeable effects 
on the growth rate of the error). It would be interest- 
ing to combine the 11- and T-formulations by using the 
former in the inner spheroidal patch and the latter in 
the other patches, closely resembling the "interpolating 
formulation" introduced in [Q- 

In this work, as in |Q, we only considered fully first 
order systems. Recently, second order in space systems 
have generated more interest, both at the continuum and 
discrete level [H IH EE El IM El E IH . We intend to 
apply the overlapping grid method to such systems and, 
most importantly, to the fully nonlinear dynamical case, 
in which non trivial issues arise. Unlike in the toy model 
cases that we have investigated so far, in the nonlinear 
case one has to face the problem of tracking a suitable 
outflow surface containing the black hole singularity and 
the dynamical generation of a grid adapted to this sur- 
face. We arc currently working on a number of related 
points. 

Overlapping grids represent possibly the simplest and 
most flexible technique capable of accurately handling 
smooth and time dependent boundaries within finite dif- 
ferencing. It is our hope that it will become a useful tool 
for the the binary black hole problem. 



When no boundaries are present (or if the gridfunctions 
are periodic), we have 



{u,Dv)h = -{Du,v)h , 



(Al) 



where {u,v)h = h^-UjVj. In |2l| it is shown that, by 
appropriately modifying the difference operator and the 
scalar product near and at the boundary, one can recover 
the summation by parts rule [4^ 



{u,Dv)h = -{Du,v)h + UjVj\ 



N 
J 10 ■ 



(A2) 



We are interested in a globally fourth order accurate 
scheme. For this purpose one could use a sixth order ac- 
curate operator at the interior with a third order accurate 
modification near the boundaries, or a fourth order accu- 
rate operator in the interior with a third order accurate 
modification near the boundaries. Whereas the first op- 
erator satisfies the summation by parts rule with respect 
to a diagonal scalar product, the second one requires a 
non-diagonal scalar product 



N 



ij=0 



(A3) 



Since the global accuracy is the same, in this work we 
choose the second option as it involves a smaller stencil 
and therefore requires fewer computational resources. Its 
modification near the boundary is given by 



1 ^ 



ijUj i = 0, ...,4 



(A4) 

^JUN-J i = iV-4, ...,il(A5) 



j=o 



where 
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APPENDIX A: MODIFICATION OF THE 4TH 
ORDER ACCURATE OPERATOR 

The fourth order accurate centered difference operator 
is given by 

DUi = (-Wi+2 + 8Ui+i - 8Mi_i + Uj_2)/(12/l) 



doo = -11/6 

dai = 3 

do2 = -3/2 

do3 = 1/3 

do4 = 

d05 = 

doe = 

dio = -24(-779042810827742869-h 

104535124033147V26116897) //i 
dn = -l/6(-176530817412806109689 + 

29768274816875927V26116897) //i 
di2 = 343(-171079116122226871-h 

27975630462649V26116897) //i 
di3 = -3/2(-7475554291248533227-H 

16484642 1 8793925 V26T16897) //i 
di4 = 1/3(^2383792768180030915 
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di5 
die 

d20 

d2i 

d22 
d23 

d2i 

d25 
d26 
d30 

dzi 

d32 

das 
dM 
dss 
dae 

dAO 

di2 
dis 
da 
di5 

^46 



1179620587812973\/26116897)//i 
-1232(-115724529581315 + 
37280576429 V261 16897) //i 




= -12(-380966843+ 86315V26TT6897)//2 
^ 1/3(5024933015 + 2010631 V26116897)//2 
= -231/2(-431968921 + 86711\/26TT6897)//2 
= (-65931742559+ 12256337V26116897)//2 
= -l/6(-50597298167+ 9716873V26116897)//2 
= -88(-15453061 + 2911\/26TT6897)//2 
= 

= 48(-56020909845192541 + 

9790180507043V26116897) //i 
= l/6(-9918249049237586011 + 

1463702013196501V26116897) /A 
= -13(-4130451756851441723 + 

664278707201077 V26116897) //i 
= 3/2(-269371Q8467782666617 + 

5169063172799767^26116897) //i 
= -1/3(6548308508012371315 + 

3968886380989379 V261 16897) //i 
= 88(-91337851897923397 + 



19696768305507\/26116897)//i 
^ 242(-120683+ 15\/26116897)//3 

264(-120683+ 15V26116897)//3 
= l/3(-43118111 + 23357^261 16897)//3 
= -47/2(-28770085+ 2259\/26Tl6897)//3 

-3(1003619433+ 11777V26116897)//3 
= -ll/6(-384168269+65747V26116897)//3 
= 22(87290207+ 10221\/26116897)//3 
= -66(3692405 + 419\/26Tl6897)//3 



APPENDIX B: HIGHER ORDER ACCURATE 
DISCRETIZATIONS NEAR THE AXIS OF 
SYMMETRY 

To illustrate the difficulties that arise with higher order 
discretizations of axisymmetric systems we consider the 
polar wave equation written in first order form 



dtP 



dpT 



(Bl) 
(B2) 



where P\p=o — 0. This system admits the following con- 
served energy 



E = 



(T^ + P'^)pdp 



(B3) 



in the sense that its time derivative gives only boundary 
contributions 

^ = 2(pTP)|p_. (B4) 

From the vanishing of P at p = we have that 

lim -dp(pP) = 2dpP (B5) 

Hence we consider the following semidiscrete approxima- 
tion 



= inn 



z = l,...,7V 
l,...,iV 



and discrete energy 

N 



(B6) 
(B7) 
(B8) 

(B9) 



and 



/i = 

/2 = 



-56764003702447356523 
+8154993476273221V26116897 
-55804550303 + 9650225V26116897 
3262210757+ 271861 V26116897 



The expression for the cr.y coefficients can be found in 
|2l| , where stability proofs for linear hyperbolic problems 
without corners can be found. 



Let us assume that Z? is a finite difference operator sat- 
isfying the summation by parts rule ljA2|l with respect to 
a diagonal scalar product. The time derivative of ljB9|) 
gives 



N 



-E = 2J2iT^D{pP), + P,p,DT,)a^h + AaToDPoh'' 



2TnPnPn + 2ro {2aDPoh - <7oD{pP)o) h . 



In order to have an energy estimate consistent with the 
one of the continuum we need to ensure that the last term 
vanishes, i.e., the difference operator at the axis has to 
satisfy 



2aDPoh = aaD{pP)o 



(BIO) 
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for some a. As was pointed out in [ij, if a first order 
accurate one sided difference operator is used on the axis, 
then IjElOp is satisfied for a = o-q/S = 1/4. 

Using Maple we were able to obtain a (non unique) 
locally second order accurate modification of the differ- 
ence operator which leads to discrete energy conserva- 
tion. However, we had to exploit the regularity condi- 
tions at the axis (P is odd and T and pP are even func- 
tions of p). The semidiscrete system 



ergy 



d 

T 

dt ' 



dt^"" 



di^' 

d-t' 

^P 
d-t' 

d 

d 

dt""^-^ 
^P 



(3Pi+4P2-3P3)//i (Bll) 

-((pP)2 + 2(pP)3)/(ll/j) 
Pi 

-(-3(pP)i + ll(pP)3 - 2{pPU)/{16h) 

P2 

-(-6(pP)i-ll(pP)2 
P3 

+ 16(pP)4 - 2(pP)5)/(26/i) 



-D{pP), 

Pt 



4,...,7V-4 



(8(pP)Ar_5 - 64(pP)Ar_4 + 59(pP)7V-2 

~3(pP)^)/(98/i) 

{%{pP)N-A - 59(pP)Ar_3 + 59(pP)7V-l 
-8(pP)Ar)/(86/l) 

{{pP)N - {pP)N-2)/{2h) 

{3{pP)N-3 + 8{pP)N-2 - 59(pP)jV-l 

+A8{pP)N)/{Mh) 
i-3To + T2 + 2T3)/{nh) 

(-GTo - 3ri + llTg - 2T4)/(16/i) 

(3ro - STi - IIT2 + 16^4 - 2r5)/(26/i) 

DTi i = 4, . . . , iV - 4 

(8rAr_5 - 64TJV-4 + 59rAr_2 - 3Tjv)/(98/l) 
(8rAr_4 - 59TAr_3 + 59Tn-1 - 8Tjv)/(86/l) 

[Tn - TN-2)/{2h) 

{3Tn-3 + 8Tn-2 - 59Tn-i + A8TN)/{Mh) 

(B12) 



where D is the standard fourth order accurate centered 
approximation of dp, conserves the following discrete en- 



E 



N 

E 

i=l 



{Tf + Pf)p,a^h + T^aoh^ 



where a = { i , ■ 
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-E = 2TnPnPn . 
dt 



(B13) 



Wc have 



(B14) 



Numerical experiments done at high resolutions suggest 
that the overall order of accuracy of the scheme is three. 

We now compare this discretization with a fourth order 
accurate approximation, obtained by simply exploiting 
the regularity conditions of the fields. We use the ap- 
proximation (|B6HB8p where D is the standard centered 
fourth order accurate finite difference operator, appro- 
priately modified at the outer boundary and use the fact 
that T and pP are symmetric across the axis. We see 
that the time derivative of the discrete energy ljB9|l is 

d 
dt 



2TmPnPn - ^(SToPi ~ 2T0P2 - 2TiPi)h 



4a 



8Pi)h 



Clearly, there is no value of a that leads to cancellations 
of the terms near the axis. Furthermore, if corrections 
of the form T^h'^ or Pl^h'^ are added to the energy, prod- 
ucts like TiPi+2 or PiTi+2 appear in the estimate. Also, 
any mixed product PiTjh? would give rise to P^ and 
terms. This discretization does not conserve the energy 
(|B9(I . In particular for a = 1/4 one is left with the axis 
terms 



(ToP2 + 2riPi)/i. 



(B15) 



Clearly, the fourth order accurate discretization based 
on regularity conditions, although not energy conserving, 
is much simpler. The introduction of ghost zones makes it 
numerical implementation trivial. Fig. [^illustrates how, 
for a particular problem and with some particular initial 
data, the energy of the fourth order accurate operator os- 
cillates about the value of the energy conserving scheme. 

Furthermore, no long term growth of the error was 
observed with the simpler fourth order approximation. 
In fact, our experiments indicate that, in general, the 
errors are much smaller than with the energy conserving 
scheme (|B11() . This is due to the lower order of local 
accuracy near the axis of symmetry and is illustrated in 
Fig.Uni 

In our overlapping grid code we adopted the simpler 
fourth order accurate approximation near the axis. 



APPENDIX C: A COMPARISON BETWEEN 
TWO FIRST ORDER REDUCTIONS 

In we performed long term stability tests to check 
that the interpolation process did not introduce un- 
wanted growth of the error. There we used M = a = 
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(3 = and had only the spherical inner patch and the 
cylindrical patch. Here we repeat the same test, but 
with nontrivial values for the black hole mass and spin. 
Interestingly, this test revealed that the approximation 
for the initial-boundary value problem for the reduction 
0-®, which we will refer to as T-reduction from now 
on, suffers from exponential growth of the error when- 
ever the domain contains a region in which the system is 
only strongly hyperbolic. These errors are concentrated 
near the excision region and, at the resolutions that we 
typically use, become evident only after about a hundred 
crossing times. 

To exclude the possibility of the overlapping grid 
method being responsible for this growth, we eliminated 
all grids, except the inner one. The growth was still 
there. Furthermore, the growth was present also when 
the second order accurate, energy conserving method for 
axisymmetric systems developed in |l4j and the modifi- 
cation of the dissipative operator constructed in were 
used. 

We recall that the first order formulation {TJ-lPl) is 
symmetric hyperbolic only outside the ergorcgion (out- 
side the black hole in the non rotating case). It is impor- 
tant to realize that the use of the discrete energy method 
to obtain stable discretization relies on the fact that the 
energy is a positive definite quadratic form of the main 
variables. Not having symmetric hyperbolicity on the en- 
tire domain, we have no guarantee of obtaining a stable 
discretization. 

This prompted us to investigate a different first order 
reduction, which, unlike the T-reduction, is symmetric 
hyperbolic everywhere. The energy associated with the 
symmetrizer, however, is not conserved and therefore we 
have no preferred way of discretizing the system. After 
describing this alternative reduction we experimentally 
compare the stability properties of the two formulations. 

As in Sec. IIIII wc start with the wave equation around 
a curved background 



. -9 ' d<i> 
Using the ADM decomposition of the metric 

ds"^ = g^^dx^dx" = -a^df + h,j{dx' + I3'dt){dx^ + dt) 

(CI) 

the wave equation becomes 

where h = det(/iy). 

We introduce the auxiliary variables 11 = {f3^di^ — 
dt^)/a and di = 9;$. Defining 



aK 



1 

7^ 



Vh 

the wave equation can be written in the form 



dt 
dU 

'dt 



ddj 
dt 



= ~ ah'^djdi - aK'^T^^jdk 

dV 



-h'^d,adj+aKU 



(C2) 



(C3) 



f3^djd, - a^^n + 8^(3^ dj - d,aT\ (C4) 



We will refer to this first order reduction of the wave 
equation as the Il-rcduction. In the particular case of 
the Kerr-Schild metric we have 



H = 



Mr 



= (l + 2iJ) 

= 



-1/2 



a' 

Vh EE det(/iy )^/^ = a^VsL sinf 



= 



1 



The principal part is given by 



/?" 
^" -a/i^" 
-an, 



where /?" = jS'^rii and /i^" = W^rii. The characteristic 
variables and speeds in the direction n are 



w 



(0;n) 
„(0;n) 



■dA : 



(3" 



where A is a vector orthogonal to n. The inverse trans- 
formation is 



n = 




dn = 


1 




dA = 


{0:n 

w^' ■ 


$ = 


{0:n 
Wis, 



(+;«) 



w 



(-;«) 



The most general symmetrizer is 

^00 
H{t, f ) = ( ?7 
rjh'^ 
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where = ^{t, x) > and rj = r]{t,x) > 0. This reduc- 
tion, unhke the T-reduction, has an everywhere positive 
definite symmetrizer (it is symmetric hyperboHc) if and 
only if the t = const, slices of spacetime are space-like. 

Using the regularity conditions (the fields $, 11, dr are 
even and dg is odd) and the fact that on the axis dga ^ 
and ah^^T^jdg = —ah^^dgde we obtain the following 
regularized equation on the axis 

dtU = (3''drU-ah''''drdr-2ah'^''dgde+ah'^rijdr 
— h^^drCedr + aKIl + 

We compare single grid, long term stability of the T- 
and Il-formulations using the same second order accu- 
rate finite difference operator modified as in 0|, where 
one sided difference operators were used at the physi- 
cal boundaries and at the axis. The dissipative oper- 
ator is modified according to prescription given in [4^ 
at the physical boundaries and regularity is used at the 
axis. No special grouping of variables is done for the 
If-formulation. Long term stability tests were done us- 
ing the forcing solution of Sec. IVllI with n = 2. At the 
outer boundary r = rmax we give data to the incoming 
characteristic variables of (In the T- formulation we 
actually give data to the incoming characteristic variables 
of HA^. One can show that both methods control the 
boundary term arising in the energy estimate |46l|.') 

Our tests revealed that the Il-formulation, despite not 
being energy conserving, outperformes the T-formulation 
in terms of stability and long term stability. Fig. 
clearly shows that if the excision boundary is placed well 
inside the event horizon the T-discretization is unstable, 
whereas the Il-discrctization is stable. Interestingly, in 
the highly spinning case, the rate of growth of the er- 
ror in the T-formulation does not seem to increase with 
resolution. Note that by excising at rmin = M = \, 
the region of lack of symmetric hyperbolicity is actually 
larger in the non-spinning case. 

In an attempt to understand the different behavior of 
the T- and Il-rcductions of the wave equation, we per- 
form a Laplace-Fourier analysis of an analogous constant 
coefficient toy model initial-boundary value problem. 

Consider the wave equation around Minkowski in 
Cartesian coordinates 

^^^'''''^,,^,,<i> = (C5) 

where rj^^ ^ — diag{ — 1, +1, +1, +1}. Under the following 
change of coordinates 

t = t' , x' = x" - (3H' (C6) 

where are constant, the wave equation takes the form 

= 2P'dtd,(l) + {S'^ - f3'p^)d,djc^ (C7) 

We consider two different types of reductions: the T- 
reduction 

dt<f> = T (C8) 

dtT = 2(5' d,T + {5'^ - l3'-p^)ddj (C9) 

dtdj = djT (CIO) 



and the Il-reduction 

dtcj) = p'd^c^ + n (cii) 

dtli = l3'diR + 6'^ddj (C12) 
dtdj = j3'ddj+djli (C13) 

Notice that the constraints Cj = dj — djcj) = 0, which are 
introduced in the reduction process, propagate differently 
in the two formulations. Whereas in the first case we have 
dtCj = 0, in the second case we have dtCj = (i'^diCj. In 
both cases the cj) variable decouples from the system. We 
will drop it in the analysis that follows. 

The T-reduction has the following symmetrizer 

+ (5''^ - I3'(3^)didj (C14) 

and characteristic speeds 

0, ^" ± 1 . (C15) 

The Il-rcduction has a simpler symmetrizer 

+ 5'^d^dJ (C16) 

and the characteristic speeds are given by 

/?",/?" ±1. (C17) 

Most importantly, whereas (|C16|I is positive definite for 
any HC14p is positive definite if and only if 5^/3* < 
1. The last condition is equivalent to the requirement 
that the vector field dt be time-like. 

Assume /3i > 1 and consider the quarter space problem 
X > with periodic solutions in y and z and no bound- 
ary conditions at a; = 0. With the Il-reduction energy 
estimates can be obtained in the outflow case using the 
energy method. For the T-formulation a Laplace-Fourier 
analysis yields the following eigenvalue problem (a sys- 
tem of ordinary differential equations in the variable x) 

dxT = sdi 
{/3f-l)d,di = 2/3i(s-z/3^c^^)di 

If an eigenvalue with 5R(s) > exists, then the initial- 
boundary value problem is ill-posed in any sense |2lj |. 
The eigenvalues are 

A± = /3i(s - iPauj'') ± ^J{s - iPauj^Y - {PI - l)wAU^^ 

(C18) 

Because of the requirement that the solution belongs to 
^2(0, -(-00), we must discard those eigenvalues which have 
positive real part. Since 5R(A±) > if 5R(s) > 0, the 
problem is not obviously ill-posed. 

This analysis reinforces our suspicion that the observed 
exponential frequency dependent growth is merely due to 
the discretization. More work is needed to exactly estab- 
lish the cause of the growth. It is possible that applying 
the Laplace transform method to the semi-discrete prob- 
lem may shed some light. 
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FIG. 7: The convergence rates Qh on the three grids obtained 
using the forcing solution w{t,r,0) with n = 2. The domains 
and coarsest resolutions are: [0, 10] x [—10, 10], 160 x 320 for 
the main grid; [1,3] x [0,7r], 80 x 240 for the inner spherical 
grid; [8, 14] x [0, tt], 240 x 240 for the outer spherical grid. The 
Courant factor is set to 0.4 and the dissipation parameter to 
a = 0.0075. The values of other parameters are: M = 1, 
P = -0.85, a = 0.99. 
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FIG. 8: This figure illustrates the importance of artificial dis- 
sipation when using overlapping grids. The discrete 1/2 norm 
of the error on the main grid at three different resolution 
for M = a = /3 = 0. For this test second order accurate 
operators, satisfying summation by parts on each grid, were 
used, but no artificial dissipation was added. The Courant 
factor is set to 1. The domains and coarsest resolutions are: 
[0,10] X [-10,10], 80 X 160 for the main grid; [1,3] x [0,7r], 
40 X 120 for the inner spherical grid; [8, 14] x [0, tt], 120 x 120 
for the outer spherical grid. The lack of stability is evident. 
On the other hand, as shown in the inset, when artificial dis- 
sipation is added to the overlapping grids scheme (cr = 0.02) 
stability is restored. 
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FIG. 9: We compare the energy conserving discretization 
IBIH with the simpler fourth order accurate one based on 
regularity. We plot the discrete energy <B9ll with q = 1/4 
as a function of time. The grid consist of 200 gridpoints and 
the domain is < p < 20. We use 4th order Runge-Kutta 
with a small Courant factor, At/ Ap = 1/10 to be close to the 
semidiscrete approximation. We use initial data T = sin®(p) 
for < p < TT and zero elsewhere. The energy of the simpler 
approximation oscillates about the value of the energy of the 
energy conserving scheme. 

0.03 r 




FIG. 10: We compare the numerical errors of the energy 
conserving discretization with the fourth order accurate one. 
We use 50 gridpoints and the domain < p < 5. The 
Courant factor is At/Ap — 1. We use the forcing solution 
method with exact solution <I> = p sin p sin t and forcing term 
F = — sin t (sin p+3p cos p) / p. The errors of the simpler fourth 
order accurate approximation are clearly much smaller than 
those of the energy conserving scheme. 
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FIG. 11: Long term stability test comparing two different re- 
ductions of the wave equation (T and 11) at different resolu- 
tions using the second order accurate discretization. We used 
the following parameters: rmin ~ 1, ^max ~ 10, tmax = 5 ■ 10'*, 
k/h = 1.5, a = 0.02, Af = 1, a = (top) or a = 0.99 (bot- 
tom), and n — 2. 
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